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RESEARCH  ARTICLE 

High-Order  CESE  Methods  for  Solving  Hyperbolic  PDEs 

David  L.  Bilyeu*a’6,  Yung- Yu  Chena,  S.-T.  John  Yu*a  and  Jean-Luc  Carnbier6 
aThe  Ohio  State  University,  Columbus,  OH  43210,  USA;  b Air  Force  Research  Lab, 

Edwards  AFB,  CA  93524,  USA 

( Received  00  Month  200x;  final  version  received  00  Month  200x) 

In  the  present  paper,  we  extend  Chang’s  (Chang  (2010))high-order  method 
for  system  of  linear  and  non-linear  hyperbolic  partial  differential  equations.  A 
general  formulation  is  presented  for  solving  the  coupled  equations  with  arbi¬ 
trarily  high-order  accuracy.  To  demonstrate  the  formulation,  several  linear  and 
non-linear  cases  are  reported.  First,  we  solve  a  convection  equation  with  source 
term  and  the  linear  acoustics  equations.  We  then  solve  the  Euler  equations  for 
acoustic  waves,  a  blast  wave,  and  Shu  and  Osher’s  test  case  for  acoustic  waves 
interacting  with  a  shock.  Numerical  results  show  higher-order  convergence  by 
continuous  mesh  refinement.  The  new  high-order  CESE  method  shares  many 
favourable  attributes  of  the  original  second-order  CESE  method,  including: 

(i)  compact  mesh  stencil  involving  only  the  immediate  mesh  nodes  surround¬ 
ing  the  node  where  the  solution  is  sought,  (ii)  the  CFL  stability  constraint 
remains  to  be  the  same,  i.e.,  <  1,  as  compared  to  the  original  second-order 
method,  and  (iii)  shock  capturing  capability  without  using  an  approximate 
Riemann  solver. 

Keywords:  CESE;  Higher  Order;  Arbitrary  Order;  CFD;  Euler  Equation;  Acoustic 
Equation 


1.  Introduction 

In  this  work,  we  extend  Chang’s  fourth-order  CESE  method  (Chang  (2010))  for  a  single 
non-linear  hyperbolic  equation  to  a  system  of  coupled  hyperbolic  partial  differential  equa¬ 
tions  (PDEs).  The  new  formulation  is  general  and  can  be  used  to  achieve  an  arbitrarily 
order  of  convergence.  To  demonstrate  the  capabilities  of  the  new  scheme,  we  apply  the 
method  to  solve  three  sets  of  equations:  (i)  the  one-dimensional  Euler  equations,  (ii)  the 
linearised  acoustic  equations,  and  (iii)  a  convection  equation  with  a  source  term. 

The  original  second-order  CESE  method  of  (Chang  (1995))  solves  the  hyperbolic  PDEs 
by  discretizing  the  space-time  domain  by  using  the  conservation  elements  (CEs)  and  so¬ 
lution  elements  (SEs).  The  profiles  of  unknowns  are  prescribed  by  assumed  discretization 
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inside  SEs.  Aided  by  the  approximation  for  the  unknowns  in  the  SEs,  space-time  flux 
conservation  can  be  enforced  over  each  CE.  The  calculation  of  space-time  flux  conserva¬ 
tion  results  in  the  formulation  for  updating  the  unknowns  in  the  time  marching  process. 
The  special  features  of  the  CESE  method  include:  (i)  The  a  scheme,  the  core  scheme 
of  the  CESE  method,  is  non-dissipative.  (ii)  The  CESE  method  has  the  most  compact 
mesh  stencil  possible,  involving  only  the  immediate  neighbouring  mesh  points  that  sur¬ 
round  the  node  where  the  solution  is  sought,  (iii)  The  method  uses  explicit  integration 
in  time  marching.  The  stability  criterion  is  CFL  <  1.  (iv)  No  Riemann  solver  is  used 
and  the  scheme  is  simple  and  efficient.  The  CESE  scheme  has  been  previous  used  to 
solve  a  multitude  of  different  physics  in  one,  two  and  three  dimensions,  including:  MHD 
(Zhang  et  al.  (2006)),  Navier-Stokes  (Zhang  et  al.  (2000),  Chang  (2007),  Venkatachari 
et  al.  (2008)),  waves  in  solids,  (Chen  et  al.  (2011),  Yang  et  al.  (2011)). 

This  paper  is  organized  as  the  following.  Section  2  reports  the  fourth-order  CESE 
method  for  the  coupled  equations  formulated  in  a  vector  form.  Section  3.3  shows  the 
application  of  the  general  formulation  to  the  one-dimensional  Euler,  linear  acoustic,  and 
convection  equation.  Section  3  provides  the  results  and  discussions  and  our  conclusion 
in  section  4. 


2.  Arbitrary-Order,  One-Dimensional  CESE  Method 

Consider  a  system  of  coupled  convection  equations: 

dU  dF 

dt  dx 


(1) 


where 


U  (ui,u2,u3r  • 

•  •  i  Um)T, 

-  conserved  variables 

Fd=  (/i,/2,/3,-- 

■  Jm)T, 

-  fluxes 

c  def  , 

»  =  («1,  S 2,  S3,  •  • 

■  ,sm)T, 

-  source  terms 

There  are  M  equations  to  be  solved  in  the  system  Eq.  (1). 

The  space-time  stencil  used  in  this  derivation  is  the  same  as  that  reported  by  Chang 
(Chang  and  Wang  (2003))  and  is  repeated  here  for  completeness.  In  Fig.  (1)  the  solid  dots 

A,  C,  and  E  are  the  solution  points.  A  is  at  ( Xj,tn )  and  C  and  E  are  at  (xj_  1/2,  tn_1^2) 
and  (xj+i/2,  tn-1/2),  respectively.  P+  is  between  M+  and  F.  P~  is  between  M~  and 

B.  The  distance  between  P±  and  M±  is  determined  by  a  parameter  r.  The  rectangles 
ABCD  and  ADEF  are  basic  CEs  (BCEs),  while  the  rectangle  BCEF  is  the  compound 
CE  (CCE)  associated  with  the  solution  point  A. 
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B  P~  M~  A  M+  P+  F 


To  facilitate  the  discussion,  we  let  SE (j,  n)  denotes  the  SE  located  at  x  =  Xj  and 
t  =  tn.  To  denote  high-order  derivatives,  we  use  the  following  notations: 

_  da+hum 
UmxHb  dxadtb 

In  SE(j, n),  the  unknown  variables  urn ,  m  =  1  are  approximated  by  a  Taylor 

series: 


u. 


i(x,t]j,n )  =f  Y 


Nm  Nm—u 

E 

a=0  6=0 


(u 


mxatb ) ; 

~Abl 


(x  —  Xj)a  ( t  —  tn) 


(2) 


where  Nm  is  the  desired  order  subtracted  by  1,  e.g.,  for  the  fourth-order  scheme,  Nm  =  3. 
The  superscript  *  represents  the  numerical  approximation  of  the  variable.  Inside  of  a  SE 
umxatb  are  constant.  The  flux  functions  /m,  m  =  1, . . . ,  M,  can  also  be  represented  with 
the  Taylor  expansion  as: 


Nm  Nm—u 


fm(x,t;j,n)  =  Y  Y 


a 


a= 0  6=0 


(3) 


Inside  a  SE,  frnx«t<‘  are  constant.  An  advantage  of  a  Taylor  series  is  that  its  derivatives 
can  also  be  expressed  as  a  Taylor  series 


A  B—a 


Ajmxzt 


(x,t',j,n)  =  Y  Y 


(Nr 


a  +  z-j-b+i 


a= 0  6=0 


albl 


(x  —  Xj)a  (■ t  —  tn) 


b  .  A=Nm~Z 
’  B=Nm—z—i 


(4) 


and 


I  mxzt 


A  B—a 

(x,  t;j,  n)  =  Y  Y 

a=0  6=0 


(/, 


mxa+ztb+‘l 


alb! 


(x  —  Xj)a  (■ t  —  tn)b 


A=Nm—z 

B=NM—z—i' 


(5) 


Equations  2  and  3  contain  2  En==L+1  n  unknowns.  In  the  following  derivation,  it  will  be 
shown,  that  for  a  given  SE(j,  n),  the  only  independent  variables  are  the  spatial  derivatives 
(' Um)j  )  {umx)j  ,  .  .  .  ,  {umxNM  ,  771  =  1,  .  .  .  ,  M. 
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We  will  now  show  two  different  methods  to  derive  the  fluxes.  The  first  method,  shown 
below,  uses  the  chain  rule  and  finds  the  derivatives  of  the  fluxes  with  respect  to  the  flow 
variables.  The  second  method,  shown  in  section  3.3  is  specific  to  a  particular  flow  physics 
but  is  easily  expanded  to  higher  orders. 

Since  it  is  assumed  that  the  fluxes  are  known  functions  of  the  flow  variables,  the  flux 
terms  in  Eq.  (3)  can  be  determined  from  the  chain  rule.  To  proceed,  we  define: 

def  dfm  r  def  d  frn  def  93fm.  (n 

J  TTll  r\  5  J TTll  k  r\  5  J TTll  k  ,p  1  *  *  * 

OUl  OUlOUk  OUlOUkOUp 

where  m,  Z,  k,p  =  1, . . . ,  M.  For  SE(j,  n),  we  obtain: 


O  £  Ncq  £  O 

OJm  \  ^  OJm  OU[ 

dyi  ^  dui  dyi  Vl  X’ 

d'2fm  _  dfm  d2Ui  d2fm  dui  duk 

dyidy2  ^  dui  dyldy2  ^  duiduk  dyi  dy2  Vl'V2 


(,T,x),  (t,t),  (x,t) 


d3fm  _  dfm  d3ui  +  d2fm  f  d2ui  duk  +  d2ui  duk  +  d2ui  duk\ 

dy\ dy2dy:i  ^  dut  dyidy2dy3  ^  duiduk  \dyidy2  dy3  dy±dy3  dy2  dy2dyi  dyx  ) 

y?  d3fm  dui  duk  dup  _  (x,  x,  x),  (t,  t,  t), 

fr^  duidukdup  dyi  dy2  dy3  yi’y2’y3  (x,x,t),  (x,t,t), 

(7) 


for  m  =  1, . . . ,  M.  Equation  Eq.  (7)  shows  the  derivatives  required  by  the  fourth  order 
scheme  but  these  equations  will  continue  up  to  the  derivatives  required  by  the  desired  or¬ 
der.  A  drawback  to  this  method  is  that  with  each  additional  derivative  the  flux  equations 
becomes  more  complex. 

Next  we  will  derive  the  equations  used  to  find  the  temporal  derivatives.  These  deriva¬ 
tives  are  solved  by  substituting  Eqs.  (2)  and  (3)  into  Eq.  (1)  in  a  given  SE(j,  n),  yielding 


du*m(x,  t;  j,  n)  df^( x,  t;  j,  n) 


=  Sr 


m  =  1, . . . ,  M. 


dt  dx 

Then  by  taking  derivatives  of  Eq.  (8)  in  both  space  and  time  we  get 


(8) 


«t)]  =  is™)]  - c rmx )?,  (u*mxtyj  =  (Smx)]  -  {f*mxx)i  «,«)?  =  (%,);  - 

This  result  can  be  written  in  a  more  general  form  as 


u. 


mxzV 


=  S. 


mxztx 


/* 


(9) 


for  m  =  1,  . . . ,  M,  z  =  0, . . .  Nm,  i  =  1,  • . .  Nm  —  z.  As  shown  in  Eqs.  (7)  and  (9),  the 
only  independent  variables  are  the  spatial  derivatives  for  each  governing  equation.  As 
such,  there  are  (Nm  +  1  )M  unknowns  for  M  equations  associated  with  a  mesh  point.  For 
example,  a  fourth-order  representation  of  the  Euler  equation  would  contain  4  unknowns 
per  equation,  giving  a  total  of  12  unknowns  for  the  one-dimensional  Euler  equations. 
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To  proceed,  for  each  of  m  =  1, . . . ,  M,  the  unknowns  are  categorized  as:  (i)  even-order 
derivatives,  («„)■,  . . . ,  umx2n  and  (ii)  odd-order  derivatives,  {umx)t}',  (fWrrr)", 

■  ■  -umx2n+i  n  =  ( OD  —  l)/2.  Where  OD  is  the  accuracy  of  the  Taylor  series.  In  what 
follows,  we  introduce  an  arbitrary  order  CESE  c-r  scheme  for  a  system  of  M  PDEs. 
The  c-t  scheme  uses  space-time  integration  for  the  advancing  formula  for  the  even-order 
derivatives,  while  the  odd-order  derivatives  are  calculated  from  a  central-difference- like 
procedure. 


2.1.  Even- Order  Derivatives 

It  can  be  shown  that  by  differentiating  Eq.  (8)  twice  with  respect  to  x,  we  obtain: 


f  ~  Ji  )  ,  dfmXX(x,t]j,  n )  _  d2s*m  _ 

I  o  O  )  ^5  1  *  *  J  1 

at  ox  axA 

or  in  a  more  general  form 

du*  2z(x,t\j,n )  df^z(x,t;j,n)  d2zs. 


(10) 


+ 


,  2  =  0, 1, ... ,  (Nm  —  l)/2.  (11) 


dt  '  dx  8x2z 

Consider  the  Euclidean  space  E2  with  the  coordinates  (£i,  £2)  =  (^1 1).  Aided  by  defining: 


def 


=  (fmx2*  j,  n),u*mx 2,  (x,  t;  j,  n)Y 


z  =  1, . . . ,  (Nm  —  l)/2, 


and  the  divergence  theorem,  Eqs.  (8)  and  (11)  can  be  transformed  into  integral  equations 
as: 


z  =  0,...,(Nm-  l)/2, 


(12) 


where  S(V)  is  the  closed  boundary  of  an  arbitrary  region  V  and  ds  is  defined  in  Fig.  (2). 


Figure  2.:  Space-time  integration  over  an  arbitrary  closed  domain  V . 
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We  define: 


def  d z+ku*n  /  Ax\ z  /  A A  k 
dxztk  \  4  )  \  4  ) 


,*  def  dZ+kfm  ( Z  ( k 

JmSiztk  dxHk  y  A  j  y  4  ) 


def  dz+ksm  /  Ax\ z  /  A  A  k 
dxztk  \  4  )  \  4  / 


(13) 


(14) 


(15) 


where  Ax  =  Xj+ 1/2  —  Xj-1/2  and  At  =  tn  —  A-1.  In  order  to  write  equations  more 
compactly,  any  local  constant  enclosed  within  a  square  bracket  will  be  evaluated  at 
the  location  specified  by  the  subscript  and  superscript  written  on  the  enclosing  square 
bracket,  e.g.: 


(«■ 


)n 

mxx  )  j 


,  Ax 


[ 


+  (Umxxt) 


At 


l  Ax  At 

/U"mxx  +  Lhnxxx  2  ^ mxxt  ^ 


J  3 


Aided  by  Eqs.  (13)  and  (14),  Eq.  (12)  gives: 


Amxz  )  /  —  A  II  Sm xz  dV 4" 


Ax 

i  Nm-z 

E 


2  (fc+1)! 


,  At 

^Jjmxk+Z  "+  ^  Jmxztk 


n—1/2 


+ 


1-1/2 


(-1  )kur 


-Af 

k+z  ^  ^  Jmxztk 


N  M  —z  —  1 


^  A, 

2  (2fe  +  l)! 


rc-l/2N 

1+1/2  , 

(16) 


Equation  (16)  provides  an  explicit  formulation  for  all  even  spatial  derivatives.  As  long 
as  the  highest  even  derivative  is  calculated  first  the  last  term  on  the  RHS  will  have  already 
been  calculated.  For  example,  in  a  fourth  order  accurate  scheme  the  conserved  variables 
are  um,umx,umx2,  and  umx 3.  In  this  case  (umx 2)"  will  be  calculated  first  followed  by 
(um)j  ■  It  should  also  be  noted  that  the  *  is  absent  from  the  source  term.  This  is  because 
the  source  term  treatment  varies  when  dealing  with  different  flow  physics  and  may  not 
require  a  Taylor  series  expansion. 


2.2.  Odd- Order  Derivatives 

In  order  to  compute  the  odd  derivatives  a  central  differencing  approach  is  applied  fol¬ 
lowing  the  c-t  scheme.  There  are  two  possible  formulations  for  the  odd-order  derivatives 
(i)  the  standard  c-t  scheme  which  is  applicable  if  there  are  no  discontinuities  present 
and  (ii)  a  re-weighted  c-t  minmod  scheme  which  is  used  if  there  are  discontinuities  in 
the  flow  held. 
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In  order  to  mitigate  the  dissipation  as  the  local  CFL  number  decreases  the  central 
differencing  is  applied  at  points  P+  and  P~ ,  where  P±  are  points  located  at 


x  (p+)  =  Xj  +  (1  +  r)—  =  xj+1/2  -  (1  -  t)  —  ,  (17) 


x(P  )  =  Xj  -  (1  +  r)—  =  Xj_1/2  +  (1  -  t)— .  (18) 

Where  r  is  the  absolute  value  of  the  local  CFL  number. 

First  we  define  u*n^z{P±)  to  be  the  Taylor  series  expansion  of  (um^)n  from  ( Xj,tn )  to 
x{P±).  Then  we  can  solve  for  wmSz+i  by  subtracting  u^z{P~)  from  «^S*(P+): 


/*T  (2 k  +  1)1  umx2k+1+*  (1  +  r)2fc!  (19) 

for  z  =  0,  2, 4, ...  ,  TV/if  —  1  and  m  =  1,  2, . . .  ,  M.  Since  we  can  not  calculate  u*nf  z  ( P ±) 
we  approximate  it  by  uL, (P1*1).  Where  uL, (P^)  is  the  Taylor  series  expansion  from 
(.Tj-i-x^tn-i/s)  to  ^(P^)  respectively. 

2  1 

El  /  .2fc 

(2k  +  1  yUm*2k+1+z  9  +  T)  ’ 

When  discontinuities  are  present  in  the  flow  field  a  re-weighting  and  or  limiting  of  the 
derivatives  is  required.  To  calculate  the  odd  derivatives  we  first  apply  re- weighting  then 
check  for  smoothness.  If  the  results  are  not  found  to  be  smooth  we  then  apply  a  limiter 
such  as  the  minrnod  operator  to  limit  the  derivatives. 

To  begin  we  outline  the  re-weighting  algorithm.  It  should  be  noted  that  all  previously 
re-weighting  schemes  used  in  second  order  CESE  schemes  should  be  applicable  to  the 
higher  order  CESE  schemes.  In  this  paper  the  derivation  of  the  W2  scheme  (Chang  and 
Wang  (2003))  will  be  presented.  First  the  function  W  is  given  as 

~  |z_|«  +  |z+|°‘  (20) 

To  remain  stable  in  the  presence  of  discontinuities  a  >  1.  The  odd-order  derivatives  are 
now  defined  by: 


u'mx‘(P+)  -u'mx'jP  ) 
2d  i  r) 


u. 


rx‘(P+)~ 


<(P- 


2(1  +  r) 


(^ma;2),'  —  (f^m—)z  (^ma;-2)  T  (Wm+) z  {7lmx+z^)  > 


(21) 


where 


(um.±)z  —  W±(umx_z  ,  Umx+*  i  az)i 


(22) 
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with 


def  ^  mxz_1  (R^ ) 

=  ± - - - 

1  +  T 


U, 


mxz^f 


def  .  1  ,, 


-o?  -  « 


where  is  the  Taylor  series  expansion  from  (xj±_i/ 2,  n— 1/2)  to  (a^i-r^j  n). 

For  solvers  that  are  greater  than  4thorder  it  was  found  that  the  re-weighting  scheme 
was  not  enough  for  stability.  In  this  case  we  need  to  apply  a  limiter  such  as  minrnod 
when  the  solution  is  no  longer  smooth.  We  define  are  smoothness  check  as 


dk+1um 

\ 

4/0  dkum 

dxk+1 

k + 1  dxk 

for  k  >  1 


(23) 


where  /3  is  an  adjustable  parameter  on  the  order  of  0.1  or  0.01.  If  Eq.  (23)  is  found  to 
be  true  then  the  minmod  limiter  was  applied.  This  smoothness  check  was  applied  to 
all  derivatives,  even  and  odd,  greater  than  1.  This  test  checks  to  make  sure  that  each 
successive  term  in  the  Taylor  series  is  a  correction  on  the  previous  term. 

The  minmod  operator  is  defined  as 

minmod(a,  b)  =  max(min(l,  a/b),  0).  (24) 

So  u^k  is  equal  to 

Uxk  —  'Mmx+fcminmod(umS_fc,  umS_p)  (25) 

The  above  equations  provide  an  explicit  formulation  for  the  odd  spatial  derivatives 
when  discontinuities  are  present  in  the  flow  field. 


2.3.  Numerical  Outline 

To  summarise  the  numerical  procedure  used  is  as  follows: 

(1)  First  calculate  all  of  the  Taylor  series  coefficients  for  the  flux  and  conserved 
variables  at  (n  —  1/2,  j  ±  1/2),  points  C  and  E  in  Fig.  1.  The  flux  and  temporal 
derivative  coefficients  are  calculated  by  Eqs.  7  and  8  respectively.  It  may  also  be 
necessary  to  calculate  the  source  term  coefficients. 

(2)  Calculate  all  even  derivatives  using  Eq.  (16),  starting  with  the  highest  unknown 
even  derivative. 

(3)  Calculate  the  highest  odd  derivative  via  Eq.  (21) 

(4)  Check  for  smoothness  using  Eq.  (23)  and  apply  minmod,  Eq.  (25),  as  needed 

(5)  Repeat  steps  3  and  4  until  all  odd  derivatives  are  found 

(6)  Check  the  even  derivative  for  smoothness,  Eq.  (23),  and  apply  the  minmod, 
Eq.  (25)  limiter  as  needed.  This  is  done  for  all  even  derivatives  except  the  zero 
derivative,  i.e.  k= 0  in  Eq.  (25). 
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3.  Numerical  Results 

The  following  test  cases  show  how  the  CESE  method  improves  as  the  order  of  accuracy  of 
the  method  employed  increases.  The  test  cases  used  for  convergence  tests  are:  (i)  a  simple 
convection  equation  with  a  source  term,  (ii)  acoustic  waves  modelled  by  the  linearised 
Euler  equations,  and  (iii)  acoustic  waves  modelled  by  the  non-linear  Euler  equations.  For 
all  three  cases,  we  calculate  the  order  of  accuracy  by  using  the  following  formula: 


where  fa  is  defined  as  the  difference  between  the  analytical  and  numerical  solution  and 
A Xi  is  the  grid  spacing  at  a  given  location  i.  In  all  cases  Ax  is  constant.  The  rate  of 
convergence  is  taken  as  the  slope  of  the  best  fit  line  through  the  points  (logio(Ax), 
logw(t2)  )■ 


3.1.  Convection  Equation  with  Source  Term 

The  first  test  case  is  the  convection  equation  with  a  source  term.  In  this  problem  we 
solve 


du  du 

— — b  a—  =  ao o  costa:) 
dt  dx 

—  2i r  <  x  <  27 r;  t  >  0 


where  a  and  Sq  are  constant.  Under  the  periodic  boundary  condition  the  analytical 
solution  to  this  problem  is 


u(x,  t )  =  cos(x  —  at)  +  Sosin(x). 


For  this  test  case  we  need  to  integrate  the  source  term.  Since  the  source  therm  is  only 
dependent  on  x  we  can  generate  an  analytical  results  equal  to 


JJ  S  =  aS0^-  sin(x')  1x^/2, 

JJ  Sxx  =  -aS0  (yf )  -y  sin(x')|x2_l/2  ,  •  •  • , 

J j  Sx2 n  =  (- l)n/2aSo  for  n  =  0,1, 


Nm  —  1 
2 


For  the  convergence  tests,  a  =  So  =  1  and  the  test  time  was  set  to  2.5^,  where  l  is  the 
length  of  the  domain.  In  all  calculations,  CFL  =  0.7.  Shown  in  Fig.  (3)  and  Table  (1),  the 
actual  convergence  rate  agrees  well  with  the  order  of  accuracy  of  the  scheme  employed. 
The  computational  scaling  shows  that  by  doubling  the  order  of  the  Taylor  series  the  time 
required  to  complete  a  simulation  will  increase  by  about  2 2,2 . 
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Table  1.:  Convergence  rates  of  the  numerical  solutions  of  the  convection  equation,  and 
the  averaged,  normalized  time  for  case  with  different  order  of  accuracy. 


Desired  Order 

Actual  Order 

Normalized  Time 

2 

2.00 

1.00 

4 

4.01 

5.22 

8 

7.95 

23.65 

12 

12.12 

59.04 

16 

16.05 

115.83 

20 

20.09 

190.95 

2th  Order.  • — • —  12th  Order - ^ 

4th  Order  — 9 —  16th  Order  — *- 

8th  Order  — * —  20th  Order  — *- 


Figure  3.:  The  £2  norm  of  numerical  solutions  of  the  Convection  equation  with  source 
term.  The  symbols  represent  the  actual  calculated  data  and  the  lines  represent  the  best-fit 
curves  of  the  data. 


Another  important  aspect  to  consider  is  whether  the  higher-order  resolution  is  worth 
the  additional  computational  cost.  For  this,  we  refer  to  Fig.  (4).  Shown  in  the  figures,  it 
is  more  efficient  to  use  a  higher-order  method  rather  than  increasing  the  resolution  for  a 
linear  solver. 
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CPU  time  (seconds) 


Order  - 

1 2th  Order 

Order  - 

— 0 - 

16th  Order 

Order  ~ 

- X - 

20th  Order 

Figure  4.:  The  I2  norm  versus  the  computational  time  for  the  solution  of  the  convection 
equation. 


3.2.  Linear  Acoustic  Equation 

The  second  test  case  is  the  acoustic  wave  solved  by  solving  the  linearised  Euler  equations. 
The  conserved  and  flux  values  are 


U  =  (ui,u2)T  =  (p,v)T , 

F  =  (  PooV,  —p)  =  (  PooU2,  —m  ] 

V  Poo  )  \  Poo  ) 


where  p.  U,  and  a  are  respectively  the  density,  velocity,  speed  of  sound.  The  speed  of 
sound  is  equal  to  y/'yp/p  with  7= 1.4.  The  values  with  a  subscript  00  are  mean  values  of 
the  flow  variables.  The  first-order  derivatives  for  the  flux  functions  are: 


fll  —  0,/i2  -  Poo 


/21 


^,/22  =  0. 

Poo 


Since  all  of  the  first-order  derivatives  are  constant  the  higher  derivatives  are  zero.  This 
reduces  the  calculation  of  all  fluxes  to  a  matrix  vector  multiplication. 

Under  periodic  boundary  conditions  the  analytical  solution  for  the  linearised  acoustic 
wave  equation  is 


£p, 


2nir 


P  =  Poo  +  ~ — cos  (  — —  {x  -  aoot) 


tt  tt  f  2mr 

U  =  Uoo  +  £cos  I  —j—\x  ~  aoo t) 

r  l 

for  —  <  a;  <  — ;  f  >  0 


where  n,  l,  and  e  are  respectively  the  number  of  waves  in  the  domain,  length  of  the 
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domain,  and  an  amplification  factor.  For  this  test  case  Poo=l,  Poo=T,  £=10~2,  n=l,  and 
1=2.  The  run  time  was  equal  to  4.25^-  which  allows  the  wave  to  propagate  through  the 
domain  4.25  times.  The  CFL  number  is  constant  throughout  the  domain  and  is  equal 
to  0.75.  As  seen  in  Fig.  (5)  our  desired  order  of  convergence  closely  matches  the  actual 
order  of  convergence.  Table  (2)  shows  the  desired  order  of  convergence,  the  actual  order 
of  convergence,  and  the  normalized  time.  Figure  6  shows  that  it  is  typically  more  efficient 
to  use  a  higher  order  scheme  then  to  use  more  points.  When  the  l2  norm  is  still  relatively 
high,  O.D(10~8)  it  was  found  that  8ththrough  20thorder  schemes  had  approximately  the 
same  computational  efficiency.  The  relative  numerical  cost  was  calculated  by  taking  the 
average  simulation  time  per  cell  per  iteration  for  multiple  resolutions  and  dividing  it  by 
the  cost  of  the  2nd  order  version.  The  computational  scaling  shows  that  by  doubling  the 
order  of  the  Taylor  series  the  time  required  to  complete  a  simulation  will  increase  by 
about  22  2. 


logic  (  A  x) 

2nd  Order  — • —  12th  Order  — *- 

4th  Order  — e —  16th  Order  — *- 

8th  Order  — * —  20th  Order  — *- 


Figure  5.:  The  l2  norm  of  numerical  solutions  where  points  are  actual  calculated  data 
and  the  line  is  a  best-fit  cure  of  the  data. 
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Table  2.:  Convergence  rates  for  the  linear  acoustic  solver  and  the  average  normalized 
time  for  case. 


Desired  Order 

Actual  Order 

Normalized  Time 

2 

2.03 

1.00 

4 

4.06 

3.43 

8 

8.12 

14.47 

12 

12.21 

34.33 

16 

16.48 

67.98 

20 

20.52 

116.89 

CN  £ 


CPU  time  (seconds) 


Order  — • — 

12th  Order 

Order  — : — 

16th  Order 

Order  — * — 

20“  Order 

Figure  6.:  The  i 2  norm  versus  the  computational  time  for  the  solutions  of  the  linear 
acoustic  wave  equations. 


The  above  cases  showed  that  we  were  able  to  achieve  higher-order  convergence  for 
coupled,  linear,  wave  equations. 


3.3.  Euler  Equation 

In  this  section  we  demonstrate  higher  order  convergence  test  as  well  as  our  schemes 
ability  to  capture  discontinuities. 

To  construct  an  arbitrary  order  CESE  solver  for  the  one-dimensional  Euler  equations, 
we  plug  the  following  unknown  vector  and  flux  function  vector  into  Eq.  (1): 

U  =  (u1,u2,u3)t  =  (p,pv,p/( 7  -  1)  +pv2/2)T , 

F  =  (pv,  pv2  +  p,  ( pE  +  p)v)T 

=  («2,(7-  l)it3  +  1/2(3  -  -  1/2(7-  1)«!/'ui)T’ 

where  p  is  the  density,  v  is  the  velocity,  E  is  the  internal  energy  and  7  is  the  ratio  of 
specific  heat.  For  the  higher  order  derivatives  the  order  of  differentiation  does  not  mater, 
1-6.  fmi'k  =  fnik.i  an(l  fmi,k,p  =  fnii,p,k  =  fmk, i,p  =  fmk,p,i  =  fmPtilk  =  fmPtkti  SO  Ollly  tile  first 
such  permutation  is  expressed.  One  method  to  find  the  flux  derivatives  is  to  substitute 
the  higher  order  derivatives  into  Eq.  7. 
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The  first-order  derivatives  of  the  flux  function  fi  are: 

fu  =  0,  fu  =  1,  /la  =  0, 

and  the  second-  and  third-order  derivatives  are: 

/i,,*=0,  l,  k  =  1,2, 3, 
fh,k,P=Q,  l,k,P  =  1,2,3. 

The  first-order  derivatives  of  the  flux  function  are: 

/2i  =  ^(7-3)^,  /2a  =  (3-7)^7,  /2a  =7-1, 

^  ’I/'l  t/l 


the  second-order  derivatives  are: 

/2„  = 


"3  ~ j  /2i  2  ...a 


_  (7— 3)u2 


f9  =  ILji  =o 

Jl2,2  443  J  J22, 3  u, 

and  the  third-order  derivatives  are: 

f  _  onl(7— 3)  r  _  o  (7— 3)^2  4?  _  n 

/2i,i,i  —  ^ - 7* - ’  / 21,!, 2  —  — ^ - ^3 - >  J21At 3  —  U, 

/2l,2,2  =  /2l,2, 3  =  0,  /2l,3,3  =  0, 

/22,2,2  d '  /22,2,3  d  ‘  /23,3,3  d. 

The  first-order  derivatives  of  the  flux  function  /3  are: 

4.  _  'YU2U3  ,  (7— l)uf  4.  _  7U3  3/.,  i\«2  41  _  7«2 

/3i  - - Tf  4  - >  /32  -  —  -  2  17  -  1J^,  J3a  -  — » 

the  second-order  derivatives  are: 

-f  _  o7«2«3  o(7-1)«2  -f  _  7«3  1  o(7-1)«2  t  _  7V2 

/3i,i  -  - 4—^ - ,  /3l  2  -  +  4— -s - ,  /3l  3  -  --^r, 

=  f  32,3  =  u7  ’  /3a,a=0, 


>  /2i,3  -  0, 
h3,3  =  0, 


the  third-order  derivatives  are: 

f  _  1  0M1;(7-1)  CIU2U3  f  _  r,7V3  A  (7— 1)M2  f  _  <)7“2 

j  3l,l,l  “  uf  b  uf  ’  j  3l,l,2  “  2  uf  y  Ilf  ’  /3l,l,3  —  2  uf  ’ 

hi, 2, 2  =  6(7  -  1)^1.  / 3l,2,3  =  / 3l,3,3  =  0 

/ 32,2,2  =  -3^4,  /S2,2,3  =  0,  / 33,3,3  =  0- 

Another  method  to  find  the  flux  derivatives  is  to  use  the  product  rule.  As  a  preliminary 
we  will  show  both  the  product  and  quotient  rule  for  any  derivative.  First  let  f(x,t )  = 
g(x,  t)h(x,  t)  then  any  derivative  of  /  can  be  expressed  as 


Qn+m  j 

dxndtm 


k= 0 1=0 


dn+m  k  lg  dk+lh 
l  k)  )  dxn~kdtm~l  dxkdtl 
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or  solving  for  qx„qJL  we  can  calculate  any  derivative  of  g 


Qn+mg  _  x 

dxndtm  h 


Qn-\-m  j 

dxndtm 


ee(:)(t)£SS ,s -».»«» 

k= 0 1= 0  v  7  v  7 


Now  we  can  express  any  derivative  of  the  fluxes 

dn+mfi  _  dn+mu2 
dxndtm  ~  dxndtm 

dn+mf2  ,  ,dn+mu3  3-7 
=  (7  -  IjTr-rwTEr  + 


dxndtm 

dn+mh 

dxndtm 


dxndtm 


EE 

k= 0  1=0 


n\  m 

k)\l 


7 


/  n\  ( m\  dn+m  k  lu2  dk+lv 

f—f  f—f  \k)  \  l  J  dxn~kdtm~l  dxkdtl 

7  -  1  dn+m-k-lU2v\  dk+lv 


k= 0  1=0 
Qn-\-m—k—l 


U3 


dxn~kdtm~l 


dxn~kdtm~l  J  dxkdtl 


We  still  need  to  express  the  derivatives  of  velocity  as  functions  of  flow  variables. 

Qn+mv  Qn+m  /  ^ 


dxndtm  dxndtm  \Ui 


1 

U\ 


dn+mu2 


EE 


n  \  /  rn\  dn+m  k  lv  dk+lu\ 


dxndtm  \kj  \  l  J  d xn~kdtm~l  dxkdtl 


for  (k,  l )  ^  (0, 0) 


(26) 


(27) 


Since  the  left  hand  side  is  always  the  highest  unknown  derivative  everything  on  the 
RHS  is  known.  Equations  (26)  and  (27)  represent  one  possible  formulation  and  is  not 
guaranteed  to  be  the  most  efficient.  We  used  Eqs.  (26)  and  (27)  in  our  simulations. 

To  show  higher  order  convergence  for  a  non-linear  equation  we  solve  the  acoustic 
equations.  There  are  some  problems  when  using  the  analytical  solution  because  the  Euler 
solver  will  capture  the  non-linearities  present  in  the  flow  field  while  the  “analytical” 
solution  does  not.  This  will  lead  to  increasing  errors  in  the  analytical  solution  as  Ax 
decreases.  To  mitigate  the  error,  the  perturbation  was  reduced  to  10-6.  For  this  test 
case  Poo=l/7>  Poo=b  n=2,  and  1= 4  and  the  simulation  time  is  2.5-C  The  CFL  number 
is  almost  constant  throughout  the  domain  and  is  equivalent  to  0.8.  The  convergence 
rates  are  shown  in  table  3.  From  this  table  it  can  be  seen  that  the  computational  cost 
associated  with  increasing  the  accuracy  is  approximately  squared,  e.g.  if  the  order  of 
accuracy  is  doubled  then  the  CPU  time  will  increase  4  times.  As  seen  in  fig.  8  using  a 
higher  order  is  typically  more  efficient  than  using  more  points. 


Table  3.:  Convergence  rates  of  the  numerical  solutions  of  the  convection  equation,  and 
the  averaged,  normalized  time  for  case  with  different  order  of  accuracy. 


Desired  Order 

Actual  Order 

Normalized  Time 

2 

2.00 

1.00 

4 

3.97 

2.84 

8 

7.48 

12.40 

12 

11.96 

32.30 
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2th  Order  — * —  8th  Order  — * 

4th  Order  — ® —  12th  Order  — » 


Figure  7.:  The  £2  norm  of  the  numerical  solutions  of  the  Euler  solver  for  solving  the 
acoustic  waves. 


CPU  time  (sec) 

2nd  Order  — * —  8th  Order  — * — 

4th  Order!«5f£aK  12th  Order  — o — 

Figure  8.:  The  £2  norm  versus  the  computational  time  for  the  solution  to  the  non-linear 
Euler  solver. 

Next,  we  demonstrate  the  new  high-order  CESE  method  by  examining  numerical  res¬ 
olution  for  shocks  and  contact  discontinuities.  We  will  run  two  different  test  cases  at 
various  resolutions  and  compare  the  results  with  a  converged  solution  obtained  by  us¬ 
ing  very  fine  mesh.  The  CESE  and  the  fifth-order  space  third-order  time  monotonicity 
preserving  (MP53)  method  (Suresh  and  Huynh  (1997)).  The  CESE  scheme  is  run  at 
4  different  orders,  2,  4,  6,  and  8th.  The  test  cases  chosen  are  Woodward’s  blast  wave 
problem  and  Shu  and  Osher’s  problem.  Woodward’s  blast  wave  problem  (Woodward  and 
Colella  (1984))  consists  of  two  shock  waves  of  different  strengths  heading  towards  each 
other  with  wall  boundary  conditions.  The  initial  conditions  are 
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(u,P,p) 


(0, 

1, 

103)  0  <  x  <  0.1 

(0, 

1, 

10“2)  0.1  <  x  <  0.9 

(0, 

1, 

102)  0.9  <  x  <  1.0 

0.0  <  t  <  0.038. 

The  second  test  case  is  Shu  and  Osher’s  problem  (Shu  and  Osher  (1989)),  in  which  a 
Mach  3  shock  moves  to  the  right  and  collides  with  an  entropy  disturbance  moving  to  the 
left.  The  boundary  conditions  are  non-reflective  and  the  initial  conditions  are 

,  f  (2.629369,  3.857143,  10.3333)  -1  <  x  <  -0.8 

^^-{(0,  1  +  0.2sin(57rx),  1)  -0.8  <  x  <  1.0 

0.0  <  t  <  0.47. 


For  both  simulations  a  in  Eq.  (20)  varies  by  this  equation. 


a  =  min 


\umx‘- 


+ 


Using  this  equation  the  value  of  a  approaches  zero  as  ucm^+  approaches  ucm^  _ .  (3  in 
Eq.  (23)  was  set  to  0.01. 

Shown  in  Figures  (9)  and  (10)  the  CESE  results  improve  as  the  order  of  the  solver 
increases.  As  the  order  of  the  solver  increase  it  becomes  more  difficult  to  suppress  the 
overshoot  around  discontinuities.  This  is  seen  in  the  8thorder  Shu-Osher  simulation.  It 
was  found  that  in  order  to  suppress  that  overshoot  the  rest  of  the  solution  became 
more  diffusive  then  the  second  order  CESE  solution.  When  evaluating  the  Shu-Osher 
simulations  the  MP  scheme  out  performs  the  second  order  CESE  scheme  and  compares 
favourably  to  the  4  and  6thorder  solutions.  The  8thorder  CESE  scheme  does  a  better  job 
then  the  MP  except  at  the  discontinuity  where  there  is  some  overshoot.  For  the  blastwave 
problem  very  few  differences  were  seen  between  the  different  simulations.  All  solvers  were 
able  to  resolve  the  shocks  with  one  or  two  points.  One  difference  is  that  in  the  400  point 
case  the  MP  scheme  does  a  slightly  better  job  when  x  is  greater  than  0.7  and  less  then 
0.8.  Also  in  the  800  point  case  the  4thorder  CESE  scheme  has  more  dissipation  then  the 
other  schemes.  This  could  be  caused  by  switching  to  the  minmod  limiter  too  soon. 


4.  Concluding  Remarks 

In  this  paper,  we  extended  Chang’s  fourth-order  CESE  method  for  one  convection  equa¬ 
tion  for  solving  a  system  of  coupled  hyperbolic  PDE’s  with  arbitrarily  high-order  con¬ 
vergence.  Numerical  results  show  that  the  extended  algorithm  can  achieve  higher-order 
convergence  for  both  linear  and  non-linear  hyperbolic  PDEs.  The  shock-capturing  ca¬ 
pability  of  the  new  method  was  comparable  to  that  of  the  original  second-order  CESE 
scheme  as  well  as  that  of  the  fifth-order  space  third-order  time  monotonicity  preserving 
scheme.  Further  development  of  the  high-order  CESE  method  can  benefit  from  further 
development  in  several  areas,  including  the  effect  of  different  limiters  on  the  higher-order 
derivatives,  and  the  effects  of  the  boundary  condition  treatments  and  the  source-term 
treatments  for  high-order  accuracy. 
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Figure  9.:  Plots  of  the  density  profiles  of  the  Woodwards  blast  wave  problem.  The  con¬ 
verged  simulation  was  done  by  using  the  second-order  CESE  with  20001  points.  For 
better  presentation,  only  a  subset  of  the  domain  is  shown  in  these  plots. 
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